Density-functional study of the Mott gap in the Hubbard model 
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We study the Mott insulating phase of the one-dimensional Hubbard model using a local-density 
approximation (LDA) that is based on the Bethe Ansatz (BA). Unlike conventional functionals 
the BA-LDA has an explicit derivative discontinuity. We demonstrate that as a consequence of 
this discontinuity the BA-LDA yields the correct Mott gap, independently of the strength of the 
correlations. A convenient analytical formula for the Mott gap in the thermodynamic limit is also 
derived. We find that in one-dimensional quantum systems the contribution of the discontinuity to 
the full gap is more important than that of the band-structure gap, and discuss some consequences 
this finding has for electronic-structure calculations. 
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The intricacies of the correlation-induced metal- 
insulator transition ('Mott transition') have fascinated 
and challenged physicists for many years, and a multi- 
tude of methods has been brought to bear on the prob- 
lem [1]. Here we study the Mott insulator from the point 
of view of density- functional theory (DFT), a method 
whose full potential for exploring the Mott phenomenon 
has not been widely recognized. It may be particularly 
timely to study the prospects for a DFT treatment of 
the Mott insulating phase because this is one of the 
two main phases of the one-dimensional Hubbard model 
(IDHM), the paradigmatic model for low-dimensional 
quantum systems. The recent interest in nanotechnol- 
ogy has brought quasi one-dimensional systems, such as 
quantum wires and carbon nanotubes, into the focus of 
current research in many-body physics and materials sci- 
ence. DFT, on the other hand, is the de facto standard 
method for the calculation of the electronic structure of 
materials, but its applicability to the very peculiar phases 
typical of low-dimensional quantum systems has not been 
systematically investigated. Here we report on the appli- 
' cation of a very recently developed density functional [2] 
to the Mott insulating phase of the IDHM. 

DFT is often applied in a single-particle mode, in which 
the eigenvalues of the Kohn-Sham (KS) equation are in- 
terpreted as a mean-field approximation to the quasi- 
particle energies. It is well known that this approxima- 
tion breaks down for systems with strong correlations, of 
which the Mott insulator is an example. On the other 
hand, as long as one restricts oneself to working with 
the particle density, the total ground-state energy, and 
quantities derivable from these, DFT is a rigorous many- 
body theory, whose predictive power is only limited by 



the quality of available approximations for the exchange- 
correlation (xc) functional. It is in this many-body mode 
that we employ DFT here. 

The choice of the appropriate approximation for the xc 
functional depends on the system under study. Here we 
investigate the IDHM, since this model represents a par- 
ticularly clear case of a Mott transition, free of unessen- 
tial complications: for nonzero interaction strength (U > 
0) the IDHM is a Mott insulator at half filling, whereas 
it is a Luttinger liquid (i.e., a one-dimensional metal) 
at all other fillings [3-5]. Following pioneering work of 
Gunnarsson, Schonhammer, and Noack [6,7] we have re- 
cently constructed an approximate xc functional for this 
model [2]. This functional is a local-density approxima- 
tion (LDA) that is based on the exact Bethe Ansatz (BA) 
solution of the homogeneous IDHM in the thermody- 
namic limit. For iV < 1 it takes the form 

£;f/^^^[n„[/] = 5]e<(n„C/), (1) 
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where e^^{n,U) = e^{n,U) — e<(n,0) — Un'^/i is the 
per-site xc energy of the homogeneous infinite IDHM, 
and e^(n, U) is the corresponding total energy, for which 
we have given an accurate parametrization, based on the 
Bethe Ansatz, in Ref. [2]. For > 1 a particle-hole 
transformation shows that e^f,(n, U) must be replaced by 
e^Jji, U) ~ e^f,(2 — n, U). Here we apply this functional 
to the Mott insulator. 

Our starting point is the exact expression for the en- 
ergy gap of a many-body system, A = I — A, where 
/ — E{N — 1) — E{N) is the system's ionization energy, 
A = E{N) - E{N+l) its electron affinity, and E{N) the 
ground-state energy of the A^-particle system. Hence, 
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A = E{N + 1) + E{N - 1) - 2E{N). 



(2) 



In ensemble DFT one can derive an alternative, but 
equivalent, representation for the true energy gap, 
namely [8] 
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where efc(iV) is the fc'th single-particle energy (KS eigen- 
value) of the -/V-particle system, and A^c is the disconti- 
nuity of the exchange-correlation potential [9] 
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where 5 — > 0+ is an infinitesimal change in the total par- 
ticle number A'^. According to ensemble DFT, for open 
systems the total energy for fractional particle number 
is given by a straight-line interpolation between the en- 
ergies obtained for integer particle number [9]. It is this 
fact that guarantees the equivalence of both representa- 
tions [8]. 

For the homogeneous infinite IDHM at n = 1 one can 
evaluate both Eq. (2) and Eq. (3) analytically from the 
parametrization underlying the BA-LDA. In either case 
we find 



where (i{U) is calculated from [2] 
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and Jo and Ji are zero and first order Bcsscl hmctions. 
(Here and below we take i = 1 as our unit of energy.) 

Representation (3) shows that there are in general two 
contributions to the many-body gap, one being the KS 
gap, Aks (this is the contribution obtained in a standard 
band-structure calculation), the other the correlation- 
induced Mott gap. Axe, which is a pure many-body effect. 
Unfortunately, the most commonly employed approxi- 
mate xc functional, the conventional ab initio LDA, does 
not have a discontinuity at integer N, so that A^^"^ = 



and A 
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This is the origin of the well known 



band-gap problem of the LDA [9]. The same applies 
to common gradient-corrected functionals. In the case 
of ab initio calculations one finds it thus necessary to 
reach beyond DFT and employ more or less systemat- 
ically derived correction terms, such as the scissors op- 
erator [10], LDA-hU [11], or LDA-hSIC [12]. However, 
the situation is different for the BA-LDA. Due to the 
particle-hole transformation involved in its construction 
this functional naturally has a discontinuity at half filling 
(occupation number n = 1), where the underlying ho- 
mogeneous infinite IDHM undergoes its metal- insulator 
transition. The BA-LDA is thus uniquely suited to test 



the density-functional theory of the many-body energy 
gap. 

To begin our analysis we present, in the left half of 
Fig. 1, the energy gap of a homogeneous IDHM with 
L = 10 sites, calculated with four different approaches: 
(i) BA-LDA gap calculated numerically from representa- 
tion (2) , (ii) energy gap obtained from representation (2) 
if the pseudo LDA of Refs. [6,13] (GS-LDA) is used in- 
stead of the BA-LDA (note that that functional does not 
have a discontinuity), (iii) BA-LDA gap calculated ana- 
lytically for an infinite homogeneous system from Eq. (5), 
and (iv) the exact gap, obtained by numerical diagonal- 
ization. Clearly, the continuous functional drastically un- 
derestimates the exact gap. In the right half of Fig. 1 we 
present the same results for a 150-site system. The last 
column is now missing since a system of this size cannot 
be treated anymore by numerical exact diagonalization. 
However, precisely because L is large, we now observe 
that the gap obtained numerically from the BA-LDA (the 
first bar) is close to the analytical result of Eq. (5) (third 
bar). 

The main result that follows from the data presented 
in Fig. 1 is that for U = 6 the dominating contribution 
to the many-body gap is due to the discontinuity. We 
have repeated the same analysis also for other values of 
U and L, and found that the critical value of U for which 
the discontinuity gap A^c becomes more important than 
the band-structure gap Aks decreases rapidly with in- 
creasing L, e.g. Uc{L = 6) = 5.2 and U^L = 10) = 4.2, 
while Uc{L — 500) = 1.8. For large systems and realis- 
tic values of U the discontinuity is thus the dominating 
contribution to the full gap. This observation is in agree- 
ment with that of Refs. [7] and [14], but contradicts the 
one of Ref. [13]. 

From our results it follows that the discontinuity con- 
tribution is absolutely crucial for obtaining the correct 
many-body gap from DFT. This implies that in order to 
solve the band gap problem encountered in the DFT de- 
scription of semiconductors and insulators it is not useful 
to optimize a continuous functional (e.g. by passing from 
LDA to GGA ) or to try to improve on the band- structure 
part of the calculation (which will only result in an im- 
proved value for Apcsjj rather one should search for a 
properly discontinuous functional. Interestingly, Ref. [15] 
reports quantitative agreement of the KS gap calculated 
from the B3LYP GGA functional and the experimental 
gap for a variety of semiconductors and insulators. These 
authors themselves conclude that this agreement is prob- 
ably fortuitious — an interpretation that is confirmed by 
the present analysis. 

Fig. 2 displays, for a 150-site system, the BA-LDA 
[2] gap and the pseudo-LDA [6,13] gap, both calculated 
numerically from representation (2) as a function of the 
interaction U. The two full curves are asymptotic results 
(valid in the thermodynamic limit) that can be extracted 
analytically from the Bethe Ansatz solution of the IDHM 
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in the limits U —^ oo and U —^ 0, respectively. For 

U ^oo [5,16] 

A(C/->oo) = i7-4+(81n2)/f/. (7) 

The leading-order term in this expansion, U ~ A, is also 
obtained immediately from our Eq. (5), since for U oo 
one has (3 = 1 [2]. For J7 — > 0, on the other hand [5,16], 

A(C/^0) = ^Vc7exp(^-^). (8) 

For U = Eq. (8) goes to zero. This, too, is reproduced 

by Eq. (5), since for U = one has (3 = 2 [2]. How- 
ever, the non-analyticity of the exact expression (8) is not 
properly accounted for by Eq. (5). This is not a defect 
of the LDA (which becomes exact in the thermodynamic 
limit of a homogeneous system) , but of our parametriza- 
tion of it, which does not incorporate this nonanalyticity. 
Work on improved parametrizations is in progress. 

As clearly seen in Fig. 2, the BA-LDA reproduces the 
limiting behaviours quite well, both for large and for 
small U, in spite of the intrinsic error of the imderlying 
parametrization. The pseudo LDA, on the other hand, 
is incapable of reproducing the exact gap for large U. 

In Fig. 3 we study the evolution of the various contri- 
butions to the many-body gap as an originally infinite 
homogeneous system with U = 6 is gradually made fi- 
nite. Both the many-body gap and the single-particle 
gap scale approximately linearly with inverse system size. 
At 1/L = 0, i.e., in the thermodynamic limit, the single- 
particle gap is zero. A band-structure calculation would 
predict the system to be a metal in this limit. The xc 
contribution to the gap, however, is nonzero and can be 
calculated from Eq. (5). For finite system sizes > 0) 
the boundaries give rise to quantization of the single- 
particle states and a discrete set of energy levels. The 
KS gap then becomes nonzero. The many-body gap is in 
this case calculated numerically from Eq. (2), and the xc 
gap is extracted from that result by means of Eq. (3). As 
seen in the inset of Fig. 3, the xc gap is a nearly quadratic 
function of 1/L, i.e., does not scale with the same power 
as the single-particle gap. 

Interestingly, the xc gap in the IDHM is larger than 
the single-particle gap all the way from L — > oo down 
to L ^ 0(1). Applied to realistic quasi one-dimensional 
systems, such as quantum wires and nanotubes, this ob- 
servation implies that the xc gap is of paramount im- 
portance for nanophysics and device technology. Un- 
fortunately, this contribution is completely missed by 
standard band-structure approaches, which only calcu- 
late the single-particle gap. Notable exceptions are 
the optimized-effective potential method [17-19] and the 
'electron-gas with a gap' approach of Refs. [20-22]. In 
view of our findings such methods may aquire additional 
significance for first-principles calculations of the elec- 
tronic structure of low-dimensional systems. 



In summary, we have demonstrated the applicability 
of the BA-LDA approach to the Mott insulating phase 
of the IDHM. We find that the discontinuity of the 
exchange-correlation functional makes a crucial contri- 
bution to the full gap, both for finite and for infinite 
systems, and for a wide range of values of U. A con- 
tinuous functional drastically underestimates the exact 
gap, whereas the BA-LDA yields quantitatively reliable 
results with very modest computational effort. This es- 
tablishes the BA-LDA, which was originally designed for 
the Luttinger- liquid phase [2], as a computational tool 
also in the insulating phase of the IDHM. Although the 
BA-LDA was constructed specifically for the IDHM and 
cannot be used directly with the ab initio Hamiltonian, 
we believe our findings concerning the discontinuity to 
be robust and qualitatively valid also for ab initio cal- 
culations. In particular, we expect that also for ab ini- 
tio calculations a reliable determination of the electronic 
structure of quantum wires and nanotubes can only be 
achieved if a functional with a proper discontinuity is 
employed. 

Acknowledgments This work was supported by 
FAFESF and CAPES. We thank M. F. da Silva for many 
discussions and providing us with the exact diagonaliza- 
tion program. KC thanks E. K. U. Gross and A. Schindl- 
mayr for helpful discussions. 



[1] N. F. Mott, Metal-Insulator Transitions 2nd. ed. (Taylor 
& Francis, London, 1990). F. Gobhard, The Mott Metal- 
Insulator Transition, Springer Tracts in Modern Physics 
Vol. 137 (Springer, New York, 1997). M. Imada, A. Fuji- 
mori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998). 

[2] N. A. Lima, M. F. Silva, L. N. Olivcira, and K. Capelle, 
submitted (2001) [cond-mat 0112428]. 

[3] E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 
(1968). 

[4] J. Voit, Rep. Prog. Phys. 58, 977 (1995). 

[5] P. Schlottmann, Int. J. Mod. Phys. B 11, 355 (1997). 

[6] O. Gunnarsson and K. Schonhammer, Phys. Rev. Lett. 
56, 1968 (1986). 

[7] K. Schonhammer, O. Gunnarsson, and R. M. Noack, 
Phys. Rev. B 52, 2504 (1995). 

[8] R. M. Dreizler and E. K. U. Gross, Density Functional 
Theory (Springer- Verlag, Berlin, 1990). 

[9] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, 
Phys. Rov. Lett. 49, 1691 (1982). J. P. Perdew and 
M. Levy, Phys. Rev. Lett. 51, 1884 (1983). L. J. Sham 
and M. Schliiter, Phys. Rev. Lett. 51, 1888 (1983). 
[10] G. Baraff and M. Schliiter, Phys. Rev. B 30, 3460 (1984). 
[11] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. 
Rev. B 52, 5467 (1995). V. I. Anisimov, J. Zaanen, and 
O. K. Andersen, Phys. Rev. B 44, 943 (1991). 
[12] A. Svane and O. Gunnarsson, Phys. Rev. Lett. 65 1148 



3 



(1990). Z. Szotek, W. M. Tcmmcrmann, and H. Winter, 
Phys. Rev. Lett. 72 1244 (1994). 
[13] K. Schonhammer and O. Gunnarsson, J. Phys. C 20, 
3675 (1987). 

[14] R. W. Godby, M. Schliiter, and L. J. Sham, Phys. Rev. 
Lett. 56, 2415 (1986). M. Lannoo, M. Schliiter, and 
L. J. Sham, Phys. Rev. B 32, 3890 (1985). 

[15] J. Muskat, A. Wander, and N. M. Harrison, Chem. Phys. 
Lett. 342, 397 (2001). 

[16] A. A. Ovchinnikov, Sov. Phys. JETP 30, 1160 (1970). 

[17] J. B. Kriogor, Y. Li, and G. J. lafratc, Phys. Rev. A 45, 
101 (1992); Phys. Lett. A 146, 256 (1990). 

[18] T. Kotani, Phys. Rev. Lett. 74, 2989 (1995). 

[19] T. Grabo, T. Kreibich, S. Kurth, and E. K. U. Gross, 
in Strong Coulomb Correlations in Electronic Structure 
Calculations: Beyond the Local Density Approximation, 
Ed. V. I. Anisimov (Gordon & Breach, 1999). 

[20] C. Gutle, A. Savin, J. B. Krieger, and J. Q. Chen, Int. 
J. Quantum Chem. 75, 885 (1999). 

[21] J. Key and A. Savin, Int. J. Quantum Chem. 69, 581 
(1998). 

[22] S. Kurth, J. P. Perdew, and P. Blaha, Int. J. Quantum 
Chem. 75, 889 (1999). 

FIG. 1. Many-body energy gap of the one-dimensional 
Hubbard model for [7 = 6. Left: L = 10 sites. First bar: 
calculated from the BA-LDA using Eq. (2). Second bar: cal- 
culated from the continuous pscudo LDA, again using Eq. (2). 
Third bar: calculated from the analytical BA-LDA expression 
(5), valid in the thermodynamic limit L —f oo. Fourth bar: 
exact result. Right: same as left, but with the first two bars 
now calculated for a system with L = 150 sites. As a con- 
sequence of the system size the bar with the exact result is 
missing. 

FIG. 2. Energy gap calculated from Eq. (2) versus interac- 
tion U for a homogeneous system with L = 150 sites. Squares: 
BA-LDA, circles: pseudo LDA. The full and dashed curves 
are the asymptotically valid analytical results for f7 — »• oo 
and f7 — > of Eqs. (7) and (8), respectively. 

FIG. 3. Energy gap of a finite IDHM as a function of in- 
verse system size 1/L. Upper curve (triangles): many-body 
gap A, calculated from Eq. (2). Lower curve (circles): 
Kohn-Sham gap Aks, calculated as a difference of sin- 
gle-particle eigenvalues. Middle curve (squares): xc gap A^c, 
calculated as the difference between the two other curves [i.e., 
from Eq. (3)]. Inset: xc gap on a smaller scale, showing ap- 
proximately parabolic behaviour of A^c as a function of 1/L. 
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L=10 



L=150 
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